#ifndef RSIM_QR_H
#define RSIM_QR_H

/* -------------------------------------------------------------------------- *
 * File: QR.h                                                                 *
 * Authors: Kishor Bhalerao                                                   *
 * Email : kishor8dm@gmail.com                                                *
 * Contributors:                                                              *
 * Email:                                                                     *
 * copyright (c) 2010 Authors.                                                *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify it    *
 * under the terms of the GNU General Public License as published by the Free *
 * Software Foundation, either version 3 of the License, or any later version.*
 *                                                                            *
 * This program is distributed in the hope that it will be useful, but        *
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
 * for more details. (http://www.gnu.org/licenses/)                           *
 * -------------------------------------------------------------------------- */

#include "Matrix.h"
#include "Vector.h"
#include "LinAlgDecl.h"

#include <cmath>

namespace RSIM
{
/**
@returns hypotenuse of real (non-complex) scalars a and b by avoiding underflow/overflow
using (a * sqrt( 1 + (b/a) * (b/a))), rather than sqrt(a*a + b*b).
*/
double hypot(const double &a, const double &b);

/** 

Note: This file was obtained from Template Numerical Toolkit (TNT).
Details of TNT can be found on the website http://math.nist.gov/tnt/
This file is minimally modified to make it compatible with RSIM.
The originial file name (available for download on TNT website) 
is jama_qr.h.

Classical QR Decompisition:
for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
orthogonal matrix Q and an n-by-n upper triangular matrix R so that
A = Q*R.

The QR decompostion always exists, even if the matrix does not have
full rank, so the constructor will never fail.  The primary use of the
QR decomposition is in the least squares solution of nonsquare systems
of simultaneous linear equations.  This will fail if isFullRank()
returns 0 (false).

The Q and R factors can be retrived via the getQ() and getR()
methods. Furthermore, a solve() method is provided to find the
least squares solution of Ax=b using the QR factors.  

(Adapted from JAMA, a Java Matrix Library, developed by jointly 
by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
*/

class QR {
	private:
		/** Array for internal storage of decomposition.
		@serial internal array storage.
		*/
	
		Matrix QR_;

		/** Row and column dimensions.
		@serial column dimension.
		@serial row dimension.
		*/
		int m, n;

		/** Array for internal storage of diagonal of R.
		@serial diagonal of R.
		*/
		Vector Rdiag;
	public:
		/** Create a QR factorization object for A. Matrix A is rectangular with (Rows>=Cols)*/
		QR(const Matrix &A);		

		/** Flag to denote the matrix is of full rank. 
		@return 1 if matrix is full rank, 0 otherwise. */
		int isFullRank() const;	

		/** Retreive the Householder vectors from QR factorization
		@returns lower trapezoidal matrix whose columns define the reflections
		*/
		Matrix getHouseholder(void)const;

		/** Return the upper triangular factor, R, of the QR factorization */
		Matrix getR() const;
	
		/** Generate and return the (economy-sized) orthogonal factor
		@param Q the (ecnomy-sized) orthogonal factor (Q*R=A). */
		Matrix getQ() const;

		/** Least squares solution of A*x = b
		@param B     m-length array (vector).
		@return x    n-length array (vector) that minimizes the two norm of Q*R*X-B.
		If B is non-conformant, or if QR.isFullRank() is false,
		the routine returns a null (0-length) vector.*/
		Vector solve(const Vector &b) const;

		/** Least squares solution of A*X = B
		@param B     m x k Array (must conform).
		@return X     n x k Array that minimizes the two norm of Q*R*X-B. If
		B is non-conformant, or if QR.isFullRank() is false,
		the routine returns a null (0x0) array. */
		Matrix solve(const Matrix &B) const;
}; // class QR


} // namespace RSIM

#endif

